clear
* Change directory to where the folder with the data is
cd "Data Archive TWT Predictions\KTHDP"

use DatasetKTHDP.dta,clear
sort rt, stable

*Merge Estimated CDFs
merge 1:1 subj trial using LeftChoice, nogenerate
merge 1:1 subj trial using RightChoice, nogenerate

*Calculate choice proportions
egen meanleft= mean(choiceleft), by(subject idpair)
sum meanleft
gen RTordered=rt

*Re-name variables to highlight what stimuli were displayed
 local k = 1
while `k' < 31{
 local i = 1
while `i' < 6{
 local j = 1
while `j' < 6{
capture noisily{
gen cdfAsubj`k'left`i'right`j'=cdfAsubj`k'id`i'`j'
gen cdfBsubj`k'left`i'right`j'=cdfBsubj`k'id`i'`j'
gen meanleftsubj`k'with`i'and`j'=meanleft if idleft==`i' & idright==`j'
sum meanleftsubj`k'with`i'and`j'
replace meanleftsubj`k'with`i'and`j'=r(mean)
}
local j = `j' + 1
}
local i = `i' + 1
}
local k = `k' + 1
}

*Loop over subjects and id pairs to create percentiles of the CDF related to the relative choice frequencies.
gen tmeank=.
 local m = 1
          while `m' < 31{
 local j = 1
          while `j' < 6{
 local i = 1
          while `i' < 6{

*get percentiles
capture noisily gen ttlow=RTordered if  meanleftsubj`m'with`j'and`i'>0.5 & 0.5/(meanleftsubj`m'with`j'and`i')>=cdfAsubj`m'left`j'right`i'
capture noisily	gen tthigh=RTordered if  meanleftsubj`m'with`j'and`i'>0.5 & 0.5/(meanleftsubj`m'with`j'and`i')<=cdfAsubj`m'left`j'right`i'

capture noisily replace ttlow=RTordered if  meanleftsubj`m'with`j'and`i'<0.5 & 0.5/(1- meanleftsubj`m'with`j'and`i')>=cdfBsubj`m'left`j'right`i'
capture noisily replace tthigh=RTordered if  meanleftsubj`m'with`j'and`i'<0.5 & 0.5/(1- meanleftsubj`m'with`j'and`i')<=cdfBsubj`m'left`j'right`i'

capture noisily	egen tlow=max(ttlow)
capture noisily egen thigh=min(tthigh)
		
*collapse new variables 
capture noisily gen tmeanksubj`m'left`j'right`i'=(tlow + thigh)/2
capture noisily sum tmeanksubj`m'left`j'right`i' 
capture noisily replace tmeanksubj`m'left`j'right`i'=r(mean) if subj==`m'
capture noisily	replace tmeank=tmeanksubj`m'left`j'right`i' if subj==`m'
				
*drop variables only used for the loop
capture noisily	drop ttlow tthigh tlow thigh
local i = `i' + 1
            }   
local j = `j' + 1
			}
local m = `m' + 1
}
save percentiles, replace 

use percentiles, clear
sort rt, stable
 
*Loop over the cases described in Theorem 2 (Alos-Ferrer, Fehr, and Netzer (2021), Theorem 3). Do this estimation out of sample (need to be repeated).
 local m = 1
          while `m' < 31{			
	local i = 1
          while `i' < 6{	
		local j = 1
          while `j' < 6{		  
	capture: gen pbarksubj`m'left`j'right`i'=.
    capture: gen pbarksubj`m'left`i'right`j'=. 
	local j = `j' + 1
            }  
    local i = `i' + 1
            }  
local m = `m' + 1   
}			
	
 local m = 1
          while `m' < 31{
 local i = 1
          while `i' < 6{			  
 local k = 1
          while `k' < 6{		  
		
		*case 1: py(`k')=1/2
		*if y is indifferent to the reference then x should have the same prop than with the ref.
		*gen pbark=py`i'`j' if  py`i'`k'==0.5
		*j is 5 the reference
		capture: replace pbarksubj`m'left`k'right`i'= 1-meanleftsubj`m'with1and`k' if  meanleftsubj`m'with1and`i'==0.5 

		
		*case 2: py(`k')>1/2, when they are both >1/2.
		*py`i'`j' : prop x over y;  FY`i'`j' CDF of x over z;  ref= entire dataset ordered RTs of reference (the grid); 
		
		*find the RT of y over z that corresponds to the percentile in the formula.
		*to do this use all values that are smaller than all RTs in the entire dataset and take the max that fulfills the condition 
		* then take the value of the CDF that corresponds to that value of the ordered RT and multiply it by the prop of x over z choices.  
		*gen hlpv1=py`i'`j'*FY`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with1and`k' * cdfAsubj`m'left1right`k' if tmeanksubj`m'left1right`i'>=RTordered
		egen hlpv2=max(hlpv1) 
		* predicted value for >1/2
		replace pbarksubj`m'left`k'right`i'=hlpv2 if  meanleftsubj`m'with1and`i'>0.5
		drop hlpv2 hlpv1
		}
		
		*case 3: py(`k')<1/2, when both smaller than 1/2.
	    *gen hlpv1=pn`i'`j'*FN`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with1and`k' * cdfBsubj`m'left1right`k' if tmeanksubj`m'left1right`i'>=RTordered
		egen hlpv2=max(hlpv1) 
		replace pbarksubj`m'left`k'right`i'=1-hlpv2 if  meanleftsubj`m'with1and`i'<0.5
				drop hlpv1 hlpv2
		}
					
local k = `k' + 1
}   		
local i = `i' + 1
}   
local m = `m' + 1   
}
save TWTPredictionsTh3r1, replace 
***************************************************************
***************************************************************
use percentiles, clear
sort rt, stable

 local m = 1
          while `m' < 31{			
	local i = 1
          while `i' < 6{	
		local j = 1
          while `j' < 6{		  
	capture: gen pbarksubj`m'left`j'right`i'=.
    capture: gen pbarksubj`m'left`i'right`j'=. 
	local j = `j' + 1
            }  
    local i = `i' + 1
            }  
local m = `m' + 1   
}			
	
 local m = 1
          while `m' < 31{
 local i = 1
          while `i' < 6{			  
 local k = 1
          while `k' < 6{		  
		
		*case 1: py(`k')=1/2
		*if y is indifferent to the reference then x should have the same prop than with the ref.
		*gen pbark=py`i'`j' if  py`i'`k'==0.5
		*j is 5 the reference
		capture: replace pbarksubj`m'left`k'right`i'= 1-meanleftsubj`m'with2and`k' if  meanleftsubj`m'with2and`i'==0.5 

		
		*case 2: py(`k')>1/2, when they are both >1/2.
		*py`i'`j' : prop x over y;  FY`i'`j' CDF of x over z;  ref= entire dataset ordered RTs of reference (the grid); 
		
		*find the RT of y over z that corresponds to the percentile in the formula.
		*to do this use all values that are smaller than all RTs in the entire dataset and take the max that fulfills the condition 
		* then take the value of the CDF that corresponds to that value of the ordered RT and multiply it by the prop of x over z choices.  
		*gen hlpv1=py`i'`j'*FY`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with2and`k' * cdfAsubj`m'left2right`k' if tmeanksubj`m'left2right`i'>=RTordered
		egen hlpv2=max(hlpv1) 
		***  predicted value for >1/2
		replace pbarksubj`m'left`k'right`i'=hlpv2 if  meanleftsubj`m'with2and`i'>0.5
		drop hlpv2 hlpv1
		}
		
		*case 3: py(`k')<1/2, when both smaller than 1/2.
	    *gen hlpv1=pn`i'`j'*FN`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with2and`k' * cdfBsubj`m'left2right`k' if tmeanksubj`m'left2right`i'>=RTordered
		egen hlpv2=max(hlpv1) 
		replace pbarksubj`m'left`k'right`i'=1-hlpv2 if  meanleftsubj`m'with2and`i'<0.5
				drop hlpv1 hlpv2
		}
					
local k = `k' + 1
}   		
local i = `i' + 1
}   
local m = `m' + 1   
}
save TWTPredictionsTh3r2, replace 
***************************************************************
***************************************************************
use percentiles, clear
sort rt, stable
 local m = 1
          while `m' < 31{			
	local i = 1
          while `i' < 6{	
		local j = 1
          while `j' < 6{		  
	capture: gen pbarksubj`m'left`j'right`i'=.
    capture: gen pbarksubj`m'left`i'right`j'=. 
	local j = `j' + 1
            }  
    local i = `i' + 1
            }  
local m = `m' + 1   
}			
	
 local m = 1
          while `m' < 31{
 local i = 1
          while `i' < 6{			  
 local k = 1
          while `k' < 6{		  
		
		*case 1: py(`k')=1/2
		*if y is indifferent to the reference then x should have the same prop than with the ref.
		*gen pbark=py`i'`j' if  py`i'`k'==0.5
		*j is 5 the reference
		capture: replace pbarksubj`m'left`k'right`i'= 1-meanleftsubj`m'with3and`k' if  meanleftsubj`m'with3and`i'==0.5 

		
		*case 2: py(`k')>1/2, when they are both >1/2.
		*py`i'`j' : prop x over y;  FY`i'`j' CDF of x over z;  ref= entire dataset ordered RTs of reference (the grid); 
		
		*find the RT of y over z that corresponds to the percentile in the formula.
		*to do this use all values that are smaller than all RTs in the entire dataset and take the max that fulfills the condition 
		* then take the value of the CDF that corresponds to that value of the ordered RT and multiply it by the prop of x over z choices.  
		*gen hlpv1=py`i'`j'*FY`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with3and`k' * cdfAsubj`m'left3right`k' if tmeanksubj`m'left3right`i'>=RTordered
		egen hlpv2=max(hlpv1) 
		***  predicted value for >1/2
		replace pbarksubj`m'left`k'right`i'=hlpv2 if  meanleftsubj`m'with3and`i'>0.5
		drop hlpv2 hlpv1
		}
		
		*case 3: py(`k')<1/2, when both smaller than 1/2.
	    *gen hlpv1=pn`i'`j'*FN`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with3and`k' * cdfBsubj`m'left3right`k' if tmeanksubj`m'left3right`i'>=RTordered
		egen hlpv2=max(hlpv1) 
		replace pbarksubj`m'left`k'right`i'=1-hlpv2 if  meanleftsubj`m'with3and`i'<0.5
				drop hlpv1 hlpv2
		}
					
local k = `k' + 1
}   		
local i = `i' + 1
}   
local m = `m' + 1   
}
save TWTPredictionsTh3r3, replace 
***************************************************************
***************************************************************

use percentiles, clear
sort rt, stable
 local m = 1
          while `m' < 31{			
	local i = 1
          while `i' < 6{	
		local j = 1
          while `j' < 6{		  
	capture: gen pbarksubj`m'left`j'right`i'=.
    capture: gen pbarksubj`m'left`i'right`j'=. 
	local j = `j' + 1
            }  
    local i = `i' + 1
            }  
local m = `m' + 1   
}			
	
 local m = 1
          while `m' < 31{
 local i = 1
          while `i' < 6{			  
 local k = 1
          while `k' < 6{		  
		
		*case 1: py(`k')=1/2
		*if y is indifferent to the reference then x should have the same prop than with the ref.
		*gen pbark=py`i'`j' if  py`i'`k'==0.5
		*j is 5 the reference
		capture: replace pbarksubj`m'left`i'right`k'= 1-meanleftsubj`m'with`k'and4 if  meanleftsubj`m'with`i'and4==0.5 

		
		*case 2: py(`k')>1/2, when they are both >1/2.
		*py`i'`j' : prop x over y;  FY`i'`j' CDF of x over z;  ref= entire dataset ordered RTs of reference (the grid); 
		
		*find the RT of y over z that corresponds to the percentile in the formula.
		*to do this use all values that are smaller than all RTs in the entire dataset and take the max that fulfills the condition 
		* then take the value of the CDF that corresponds to that value of the ordered RT and multiply it by the prop of x over z choices.  
		*gen hlpv1=py`i'`j'*FY`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with`k'and4 * cdfAsubj`m'left`k'right4 if tmeanksubj`m'left`i'right4>=RTordered
		egen hlpv2=max(hlpv1) 
		***  predicted value for >1/2
		replace pbarksubj`m'left`i'right`k'=hlpv2 if  meanleftsubj`m'with`i'and4>0.5
		drop hlpv2 hlpv1
		}
		
		
		*case 3: py(`k')<1/2, when both smaller than 1/2.
	    *gen hlpv1=pn`i'`j'*FN`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with`k'and4 * cdfBsubj`m'left`k'right4 if tmeanksubj`m'left`i'right4>=RTordered
		egen hlpv2=max(hlpv1) 
		replace pbarksubj`m'left`i'right`k'=1-hlpv2 if  meanleftsubj`m'with`i'and4<0.5
				drop hlpv1 hlpv2
		}
					
local k = `k' + 1
}   		
local i = `i' + 1
}   
local m = `m' + 1   
}
save TWTPredictionsTh3r4, replace 
***************************************************************
***************************************************************
use percentiles, clear
sort rt, stable
 local m = 1
          while `m' < 31{			
	local i = 1
          while `i' < 6{	
		local j = 1
          while `j' < 6{		  
	capture: gen pbarksubj`m'left`j'right`i'=.
    capture: gen pbarksubj`m'left`i'right`j'=. 
	local j = `j' + 1
            }  
    local i = `i' + 1
            }  
local m = `m' + 1   
}			
	
 local m = 1
          while `m' < 31{
 local i = 1
          while `i' < 6{			  
 local k = 1
          while `k' < 6{		  
		
		*case 1: py(`k')=1/2
		*if y is indifferent to the reference then x should have the same prop than with the ref.
		*gen pbark=py`i'`j' if  py`i'`k'==0.5
		*j is 5 the reference
		capture: replace pbarksubj`m'left`i'right`k'= 1-meanleftsubj`m'with`k'and5 if  meanleftsubj`m'with`i'and5==0.5 

		
		*case 2: py(`k')>1/2, when they are both >1/2.
		*py`i'`j' : prop x over y;  FY`i'`j' CDF of x over z;  ref= entire dataset ordered RTs of reference (the grid); 
		
		*find the RT of y over z that corresponds to the percentile in the formula.
		*to do this use all values that are smaller than all RTs in the entire dataset and take the max that fulfills the condition 
		* then take the value of the CDF that corresponds to that value of the ordered RT and multiply it by the prop of x over z choices.  
		*gen hlpv1=py`i'`j'*FY`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with`k'and5 * cdfAsubj`m'left`k'right5 if tmeanksubj`m'left`i'right5>=RTordered
		egen hlpv2=max(hlpv1) 
		***  predicted value for >1/2
		replace pbarksubj`m'left`i'right`k'=hlpv2 if  meanleftsubj`m'with`i'and5>0.5
		drop hlpv2 hlpv1
		}
		
		
		
		*case 3: py(`k')<1/2, when both smaller than 1/2.
	    *gen hlpv1=pn`i'`j'*FN`i'`j' if tmeank`i'`k'>=ref
		capture noisily{
		gen hlpv1= meanleftsubj`m'with`k'and5 * cdfBsubj`m'left`k'right5 if tmeanksubj`m'left`i'right5>=RTordered
		egen hlpv2=max(hlpv1) 
		replace pbarksubj`m'left`k'right`i'=1-hlpv2 if  meanleftsubj`m'with`i'and5<0.5
				drop hlpv1 hlpv2
		}
					
local k = `k' + 1
}   		
local i = `i' + 1
}   
local m = `m' + 1   
}
save TWTPredictionsTh3r5, replace 


clear
* Change directory to where the folder with the data is
cd "Data Archive TWT Predictions\KTHDP"

use DatasetKTHDP.dta, clear
*Loop over the 5 alternative out-of-sample estimations
 local m = 1
          while `m' < 6{
use TWTPredictionsTh3r`m', clear
*Merge estimated risk attitudes
merge m:1 subj using UtilityEstimationCARAKTHDP, nogen
egen m1 = max(out1)
egen m2 = max(out2)
gen maxOut=max(m1,m2)
*Generate choice proportions from estimated risk attitudes
gen RUMCARAleft = normal(prob1 * ( ((1-exp(-est_r_u`m'*out1))/(1-exp(-est_r_u`m'*(maxOut))))-(prob2 * (1-exp(-est_r_u`m'*out2)))/(1-exp(-est_r_u`m'*(maxOut))))/noise_u`m') if prob1<prob2
replace RUMCARAleft = normal(prob1 * ( ((1-exp(-est_r_u`m'*out1)))-(prob2 * (1-exp(-est_r_u`m'*out2))))/noise_u`m') if prob1>prob2

*Collapse the dataset at the subject and choice pair level
collapse choiceleft pbarksubj* idleft idright RUMCARAleft , by(subj idpair)

*calculate mean errors for TWT and EU
gen diff=.
 local j = 1
          while `j' < 31{	
 local i = 1
          while `i' < 6{	
 local k = 1
          while `k' < 6{
replace diff = choiceleft-pbarksubj`j'left`i'right`k' if idleft==`i' & idright==`k' & subj==`j'
local k = `k' + 1
            }   		
local i = `i' + 1
            }  
local j = `j' + 1
            }  			
gen diffRUMCARA = choiceleft-RUMCARAleft

*Calculate absolute mean error
gen distCARA`m' = abs(diffRUMCARA)
gen dist`m' = abs(diff)

collapse dist`m'  distCARA`m' , by(subj)
sum dist`m' distCARA`m' 
signrank dist`m' = distCARA`m'
save Th3CARApred`m', replace 
local m = `m' + 1   
}

*Merge the different out-of-sample predictions
use Th3CARApred1, clear
merge 1:1 subj using Th3CARApred2, nogen
merge 1:1 subj using Th3CARApred3, nogen
merge 1:1 subj using Th3CARApred4, nogen
merge 1:1 subj using Th3CARApred5, nogen

*Re-shape datasets for ease of comparisons
reshape long dist distCARA, i(subj) 
collapse dist distCARA, by(subj)
sum dist distCARA
save ChoicePredictionFechnerCARAKTHDP, replace
*Adapt this to get CRRA, CARA-RPM, and CRRA-RPM predicitions.
*saved files should be:
* ChoicePredictionFechnerCRRAKTHDP
* ChoicePredictionFechnerRPMCARAKTHDP
* ChoicePredictionFechnerRPMCRRAKTHDP
